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1.3  INTRODUCTION 


^  The  computer  program  PAPST,  standing  fot^ Plastic  Axisymmetric/Planar 
STructure^ ^Reference  ifj  has  been  developed  with  theNavy  with  the  aim  of 
/^-^addressing  the  elastic-plastic  crack  problems  which  arise  frcm  material  and 
structural  considerations  associated  with  the  use  of  high  toughness 
materials.  Both  power  hardening  and  multilinear  models  for  uniaxial 
stress-strain  response  are  available.  Enriched  elements  may  be  used  to  model 
the  material  crack  tip  response  for  the  multilinear  material  model.  For 
crack  problems,  both  path  values  for  the  J  integral  and  values  based  on  the 
amplitude  of  the  singular  solution  can  be  computed. 

s-'-  b  Xc 

The  program  has  beeh  used  in  support  of  laboratory  programs  to  develop 
material  property  data  for  J_L  and  other  pertinent  paraneters  for  fracture 
characterization.  More  recently,  the  capabilities  of  PAPST  have  been 
extended  to  include  simulation  of  quasi-static  crack  growth,  procedure 

is  of  interest  in  determination  of  the  Tearing  Modulus '(Reference  2)Nand  in 
assessing  the  stability  of  crack  growth. 


Incorporated  into  PAPST  is  its  predecessor  ApES. . stand ing  for 
Axisymmetric/Planar  Elastic^ Str uctures^tRef erences  3 , 4)7  Also  developed 
/  "with  the  Nav£7  APES^performs^l inear  elastic  fracture  mechanics  for  Mode  I 
and  Mode  II  problems.  It  has  the  capability  to  handle  crack  face  loadings 
and  a  limited  capability  to  deal  with  orthotropic  materials. 
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2.0  PROGRAM  DESCRIPTION 


This  report  is  to  review  and  update  the  theoretical  basis  and  present 
capabilities  of  the  computer  program  PAPST.  This  program,  and  its 
predecessor  APES,  have  been  in  development  for  ten  years  by  the  Navy.  The 
present  revision  of  PAPST  has  involved  a  major  restructuring  of  the  basic 
computer  code  to  facilitate  future  development.  Several  options  that  were 
in  previous  revisions  have  been  dropped  because  they  proved  to  be 
inaccurate.  Many  of  the  Important  capabilities  have  been  rederived  to  be 
internally  consistent.  New  features  have  been  added  to  enhance  the  progr^n. 
The  modifications  are  extensive,  requiring  a  complete  review  of  the  program. 

Present  capabilities  include: 

1)  Elastic  orthotropic  material  properties. 


2)  Elastic/perfectly  plastic,  power  hardening,  and  multilinear 
isotropic  models  for  elastic/plastic  material  behavior. 

3)  Cubic  isoparametric  plane  elements  with  three  (9-node)  or  four 
(12-node)  sides,  with  optional  enrichment  to  accommodate  the  crack 
tip  asymptotic  solution. 

4)  Concentrated  and  distributed  elastic  springs. 

5)  Tied  nodal  displacements. 

6)  Concentrated  and  distributed  tractions. 

7)  Crack  face  pressure  loading. 

8)  Thermal  loading. 

9)  Mode  I  and  II  elastic  stress  intensity  factors. 

10)  Mode  I  plastic  stress  intensity  factors. 

11)  J-integral  calculation  from  stress  intensity  factors  and  from  path 
integrals. 

12)  Finite  strain  treatment. 

13)  Load,  displacement  or  thermal  controlled  incremental  loading. 

14)  Wavefront  reordering. 

15)  Plotting  of  displacements,  stress  contours  and  thermal  contours. 

Both  enriched  and  regular  models  can  be  analyzed  using  the  non-linear 
stress-strain  relationships  detailed  in  Section  2.1.  The  program  will 
automatically  scale  the  elastic  solution  to  yielding  of  the  highest 
stressed  quadrature  point.  Subsequent  increments  are  either  defined  by  the 
user  or  evenly  stepped  to  full  loading.  Inclusion  of  the  large  deformation 
calculation  is  optional.  The  degree  of  accuracy  can  also  be  user  supplied. 
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2.3.1  Plasticity  Theory 


PAPST  incorporates  incremental  J_-flow  theory  plasticity  (von  Mises 
yield  criterion) .  The  components  of  the  deviatoric  strain  rate  are 
expressed  in  terms  of  the  current  deviatoric  stress  components  and  the 
components  of  the  deviatoric  stress  rate  as 


+!f  Sij  ae 


1  +  v 
E 


"°e  =  cyd] 

[fe  >  0  J 


S..  (otherwise) 
1  J 


(2.0. 1-1) 


where  e„  are  the  deviatoric  strain  rate  components 

1  *  * 

eij  eij  ”  3  "pp.  ij 


S.j  are  the  current  deviatoric  stress  components 

Sij  =  °ij  "  3  aPP 


(2. 0.1-2) 


(2.0. 1-3) 


S, .  are  the  deviatoric  stress  components  measured  from  the  center  of 
the  current  yield  surface 


s',i  ■  S,J  - 


(2.0. 1-4) 


a.,  are  the  coordinates  in  stress  space  of  the  center  of  the  yield 

surfac‘d  c  is  the  von  Mises  effective  stress 
e 


'£  c  s 
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(2. 0.1-5) 


and  a  is  the  effective  stress  as  measured  from  the  center  of  the 
yield  surffce  and  a ^  is  the  yield  stress. 
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The  function  f(oo)  in  Equation  2.0. 1-1  is  derived  from  a  unaxial 
stress-strain  curve  as  described  subsequently  (Reference  5) 


as 


The  hydrostatic  strain  rate  is  given  in  terns  of  the  mean  stress  rate 


•  -  ll hL  A 

cpp  E  pp 


(2. 3. 1-7) 


The  center  of  the  yield  surface  can  move  under  plastic  deformation  at 
a  rate  proportional  to  the  projection  of  the  stress  rate  vector  onto  the 
local  normal  to  the  yield  surface.  The  motion  of  the  yield  surface  center 
is  directed  along  the  vector  describing  the  deviatoric  stress  relative  to 
the  yield  surface  center, 

.  (2.0. 1-8) 
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10  [otherwise] 

0  <  2  <  1 


where  £  is  to  be  specified  by  the  program  user.  Isotropic  hardening  is 
achieved  by  setting  2  equal  to  zero  and  kinematic  hardening  results  when  3 
is  set  equal  to  one.  Intermediate  cases  between  these  two  extremes  may  be 
achieved.  These  possibilities  are  indicated  in  Figure  2. 0.1-1  which  shows 
various  yield  surfaces  on  the  tt  plane  in  stress  space,  and  corresponding 
examples  of  uniaxial  bilinear  stress  strain  curves  obtained  with  each 
hardening  model.  For  uniaxial  or  proportional  states  of  stress,  all 
hardening  models  would  give  the  same  result  provided  there  is  no  reversed 
yielding.  For  multi-axial  states  of  stress,  or.  the  other  hand,  the 
kinematic  and  isotropic  hardening  models  can  give  significantly  different 
results  if  the  loading  is  strongly  non-proportional  and/or  if  reverse 
yielding  occurs. 

2.3.2  Finite  Strain  Treatment 


PAPST  incorporates  finite  strain  effects  through  the  use  of  an  updated 
Lagrangian  formulation  (Reference  6),  i.e.,  the  coordinate  system  is 
convected  with  the  deformation.  In  this  coordinate  system  the  relationship 
between  true  strain  rate  and  the  deformation  rate  (or  velocity)  is 
unchanged  from  the  small  strain  theory  equations.  Also  the  true  (or 
Cauchy)  stress  rate  tensor  relates  to  the  true  strain  rate  in  the  classical 
sense , 
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ISOTROPIC  (3*0) 


KINEMATIC  (g-1) 


7r  Plane  Comparison  of  Subsequent  Yield  Surfaces  for 
Various  Hardening  Rules 


a 


Uniaxial  Bilinear  Stress-Strain  Curves  for  Kinematic 
and  Isotropic  Hardening  Models 

Figure  2.0. 1-1  Diagrams  of  Isotropic  and  Kinematic  Hardening  Models 
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where  the  subscripts  take  values  of  one  to  three  and  repeated  subscripts 
imply  summation. 


The  Cauchy  stress  tensor  is  not  invariant  under  finite  rotations,  but 
rather  at  zero  strain  rate  it  changes  at  the  rate  (Reference  31) . 


.c 
3  i  J 


=  W.  o  • 
TP  PJ 


-  W  . 

PJ  ip 


(2. 0.2-2) 


where  W. .  is  the  rotation  rate  tensor.  Equations  2. 0.2-1  and  2. 0.2-2  are 
combined-1  to  express  the  Cauchy  stress  rate  in  terms  of  the  true  strain  and 
rotation  rates  as 


,c 


=  Ci  3  kl  ^  kl 
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(2. 0.2-3) 


2.0.3  Iterative  Solution  Procedure 


The  governing  matrix  equation  for  finite  elements  is 


K  U  =  P 


(2. 0.3-1) 


Where:  K  is  the  stiffness  matrix 

U  is  the  displacement  vector 

P  is  the  load  vector 

For  elastic  solutions,  the  stiffness  matrix  is  constant,  and  the 
solution  can  be  determined  explicitly.  Incremental  plasticity  Implies  that 
the  stiffness  matrix  is  a  function  of  the  displacement  vector.  Because  of 
this,  the  stiffness  matrix  can  only  be  determined  instantaneously.  This 
requires  that  the  loading  be  incremental  and  that  the  solution  must  be 
iterated  until  the  instantaneous  stiffness  matrix  and  the  incremental 
displacements  agree  with  the  applied  loading  increnent,i .e. , 


K(U  -U  -  -P 


(2. 0.3-2) 
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The  formulation  of  the  instantaneous  stiffness  matrix  is: 


k(u) 


-  I  />  T(u)  C(u)B(u)  dA 
elements  area 


(2. 0.3-3) 


Equation  2. 0.3-2  is  solved  for  the  estimate  of  aU  by  using  the 
previous  U  vector  to  calculate  K(u) . 

The  load  correction  vector  is  then  determined  as 

2  *  -  2  KUIaU^P  (2.0 

increments  increments 


and  the  equation 


K(u)aa(J  =  aaR 


(2. 0.3-5) 


is  solved  to  obtain  the  correction  to  All . 

However,  it  is  not  necessary  to  perform  the  complete  computation 
implied  in  Equation  2. 0.3-3  to  obtain  the  load  correction  vector,  rather 


K(u)  *U  = 


2 

elements 

2 

elements 

2 

elements 


)  c  (u)  B(u)  dA  aU: 


area 


f  BT(u!  Ctu 


)  a  edA 


area 


/  BT(u 
area 


(2. 0.3-6) 


)a a  dA 


An  exception  to  this  procedure  should  be  made  for  any  elements  for 
which  the  stiffness  matrix  is  readily  available.  This  is  the  case  for  both 
concentrated  and  distributed  springs  since  they  are  elastic,  i.e., 
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K(u)aU  J  3T(u)  dA  ][  KspAU 


(2. 3. 3-7) 


elements  area 


springs 


Two  different  convergence  criteria  have  been  set-up: 


i 

d  AAPi  /i  p.|]1/2  <_  conv-j 


(2. 3. 3-8) 


for  overall  convergence,  and 


[|&AP|/|iP|  ]V2  £  conv. 


(2. 0.3-9) 


for  increment  convergence.  Both  must  be  satisfied  in  PAPST  at  each 
increment. 


2.3.4  Yield  Surface  Modification 

Cnee  an  increment  has  converged,  the  yield  surfaces  of  each  node  and 
quadrature  point  are  modified  to  reflect  the  new  stress  state.  The  new 
coordinates  of  the  center  of  the  yield  surface  are: 


,  =  aij  +  eu;  -  ri  sjj/o; 


(2. 0.4-1) 


where  r  =  radius  of  the  previous  yield  surface 


=  r  f  (1-3)  (o‘p  -  r) 


(2. 0.4-2) 


Che  yield  surface  is  modified  'when  the  increnent  calculation  has 
converged . 
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2.1  Material  Models 


Hie  material  model  describes  the  relationship  of  the  stress  vector 
to  the  strain  vector: 


a 

a. 


(2.1-1) 


£  =  C"1 0  ■  D  a 


(2.1-2) 


a  is  the  stress  vector 

a, 

^  is  the  strain  vector 

Where  C  and  D  are  the  constitutive  material  matrices 

The  material  matrices  in  the  elastic  range  are  constant  and  may  be 
specified  with  direction  dependent  (orthotropic)  properties.  In  the 
plastic  analysis,  the  matrices  are  restricted  to  isotropic  properties  but 
may  describe  a  non-linear  stress-strain  relationship.  The  material 
matrices  will  then  reflect  the  instantaneous  relationship.  Options  are 
available  for  either  power  hardening  or  multilinear  material  models. 
Either  of  these  models  can  be  used  to  approximate  a  true  stress/true 
strain  curve  that  has  been  experimentally  derived. 

2.1.1  Orthotropic  Elastic  Model 

Hie  orthotropic  stress-strain  relation  for  an  axisymmetric  element 
is: 
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The  ten  independent  variables  reduce  to  seven  if  we  impose  the 
Maxwell's  reciprocal  relations: 
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vyx  ’  vyz  vzy  *  vxz  vzx 


(2. 1.1-2) 


Isotropic  properties  are  a  special  case  in  that  there  are  only  two 
independent  parameters: 


Poisson's  Ratio  v  =  vxy  s  vyZ  =  vzx 


Young's  Modulus  E  =  E  =  E  =  Ez 
*  x  y 


(2. 1.1-3) 
(2.1. 1-4) 


G  =  Gxy  =  E/2 ( 1  +  v) 


(2. 1.1-5) 


For  plane  strain,  the  constraint  exists  that  the  out-of-plane  strain 
must  be  zero: 


e2  s  0  (2. 1.1-6) 

For  this  case  the  C  matrix  remains  unchanged,  but  the  D  matrix  becomes: 
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(2. 1.1-7) 


For  plane  stress  the  constraint  is  on  the  out-of-plane  stress: 
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(2. 1.1-8} 


The  D  matrix  remains  unchanged,  but  the  C  matrix  is  modified  similarly  to 
the  D  matrix  in  plane  strain. 
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(2. 1.1-9) 


The  program  has  fully  implemented  the  use  of  the  orthotropic 
relationship  for  elastic  problems. 

2.1.2  Elastic/Perfectly  Plastic  Model 

For  a  uniaxial  stress  state  the  following  relationship  exists  for 
this  model: 


a  - 


for 

e  '  °yd/E 

(2. 1.2-1) 

for 

E  -  °yd/E 

This  relationship  is  useful  in  plasticity  calculations  only.  It 
is  not  used  in  fracture  mechanics  or  J-integral  analyses.  The 
application  of  this  relationship  is  primarily  for  comparison  of  results 
to  theoretical  solutions  for  vdiich  it  is  analytically  convenient  to 
assune  elastic/perfectly  plastic  behavior. 

The  power  hardening  and  multilinear  models  can  be  specified  such 
that  they  simulate  elastic/perfectly  plastic  material  response. 
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2.1.3  Power  gardening  (Ramberg -Osgood)  Model 


In  the  caae  of  the  power  hardening  simulation  of  material  response 
(Ramberg -Osgoo:  >  ,  the  uniaxial  response  is  given  as 


-/£  for  o  < 


r/E  +  3[  f  (  f  ) 
E 


yd 

n-1  "yd 


-f-  1  for 


a  >  a 


yd 


(2. 1.3-1) 


The  correlatio-  oaraneters  and  the  initial  yield  stress  are  chosen  to 
best  model  the  :bserved  behavior.  Note  that  if  n  -  1,  A  bilinear 
material  model  ;s  achieved.  In  this  case  the  strain  hardening  slope  is 
given  by: 


13/m  =  E/0  +  ») 


(2. 1.3-2) 


For  the  power  -..ardening  material  model,  the  plastic  component  of  the 
strain  is  giver,  by 


ao 


ao 


& 


'P  '  Eaydn_1 


and  the  plastic  strain  rate  is  thus 


n-1 . 

ana  a 


(2. 1.3-3) 


(2. 1-3-4) 


The  power  hardening  material  model  is  illustrated  in  Figure  2.1-1.  For 
this  model,  the  function  f(  ag)  (Equation  2. 0.1-1)  which  relates  the 
plastic  strain  rate  to  plastic  stress  rate  is 


f(ae)  = 


(2. 1.3-5) 
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Power  Hardening  Material  Model 
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Multilinear  Material  Model 


Figure  2.1-1  Available  True  Stress  -  True  Strain  Models 


/ti  Arthur  IX  Little,  Inc 


ffiNUXa  iNIWMMiAA*) 


2.1.4  Multilinear  Model 


The  mathematical  form  for  the  multilinear  model  of  material  behavior 
is 


£ 


“2 

+  r  (a3 


+  (o-  om) 


(2. 1.4-1) 


a 


m 


<  a  $  o 


m+1 


(2.1. 4-2) 


I 

For  this  material  model,  the  plastic  strain  rate  is  given  by 


=  ao/E 
m  e 


(2. 1.4-3) 


and 


f(oe)  •  VEoe 


(2. 1.4-4) 


The  multilinear  material  model  divides  the  uniaxial  response  into 
one  elastic  and  m-1  plastic  straight  line  segments,  vrfiere  -m  is  used  as 
data  input  to  indicate  the  nunber  of  stress-strain  pairs  or  "breakpoints" 
needed  to  define  the  multilinear  model.  The  user  then  specifies  the  true 
stress  and  true  strain  at  yield  and  (m-1)  other  points  on  the  response 
curve.  The  last  (mth)  segment  continues  at  constant  slope  to  infinite 
strain  regardless  of  the  strain  level  actually  used  to  define  the  slope 
of  this  segment.  An  exanple  of  this  model  is  shown  in  Figure  2.1-1. 
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2. 2  Elements 


The  element  library  of  PAPST  contains  six  elements: 

1)  12-node  isoparametric  quadrilateral 

2)  12-node  enriched  quadrilateral 

3)  9-node  isoparametric  triangle 

4)  9-node  enriched  triangle 

5)  2-node  concentrated  spring 

6)  8-node  distributed  spring 

Except  for  the  concentrated  spring,  the  element  stiffness  matrix 
formulation  is  based  on  the  Gaussian  quadrature  integration  described  in 
Zienkiewicz  (Chapter  8,  Reference  7).  The  concentrated  spring  stiffness 
matrix  can  be  defined  explicitly  so  that  integration  is  not  necessary. 

The  following  discussion  is  provided  so  that  the  user  can  understand 
the  basic  theory  and  so  that  the  notation  can  be  defined. 

Isoparametric  elements  are  defined  as  elements  in  vrfiich  the  functions 
used  to  describe  the  geometry  and  the  displacements  are  the  sane.  Depending 
on  the  element,  these  shape  functions  can  describe  a  line,  an  area,  or  a 
volune.  Using  these  functions,  we  can  map  a  distorted  shape  in  physical 
coordinates  to  a  simple  shape  in  local  or  mapped  coordinates.  The  stiffness 
matrix  of  the  simple  shape  can  be  found  generically  and  used  for  a  larger 
class  of  physical  element  shapes.  For  example,  the  12-node  element  is 
applicable  to  rectangles,  trapezoids,  parallelograms,  four-sided  circular 
sectors,  and  so  on  for  any  element  with  four  sides  in  which  each  side  can  be 
described  adequately  by  a  cubic  function. 

For  purposes  of  this  discussion,  we  will  use  shape  functions  that 
define  an  area  element.  The  shape  is  restricted  to  the  X-Y  physical 
cooodinate  system.  In  mapped  space,  this  translates  to  the  s-t  local 
coordinate  system.  The  shape  functions  will  be  polyncminal  in  the  s-t 
coordinate  system. 

The  geometry  can  be  described  at  a  point  by  the  shape  functions 
evaluated  in  the  mapped  space  multiplied  by  the  nodal  values  of  the  global 
coordinates.  (Note:  there  are  as  many  shape  functions  as  nodes.) 

n  n  (2.2-1) 

X  =  ]T  Ni(s,t)X.  ;  Y  N.(s,t)Yj 
i=1  i=l 


where:  N^  are  the  shape  functions 

X.  and  Yj  are  the  nodal  global  coordinates 


/L  Arthur  D,  Little,  Inc. 


Similarly,  the  displacements  at  these  nodes  are  described: 


n  n 

U  =  2  VS*‘>U,  J  v  -  J  N1<S't)Vi  (2.2-2) 

i=l  i=l 


vtfhere  IL  and  V.  are  the  nodal  global  displacements  (U  corresponds  to  X  and  V 


i 

,  strains  within  the  area  can  be  found  by  using  the  appropriate 
derivatives  of  the  shape  functions: 


n 

«=BS  =  Ib 

i*l 


V 


i 

i . 


(2.2-3) 


where:  B  is  the  matrix  of  shape  function  derivatives  appropriate  for  the 
type  of  analysis.  In  plane  stress  or  plane  strain  the  following  holds: 


LYxy  J 


=2 

i=l 
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while  in  axisyronetric  problems: 
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Note:  in  PAPST  the  axisyrrmetric  coordinate  system  is:  r,  z,  e 
corresponding  to  (x,  y,  z) .  Also,  for  consistency  all  problems  take  the 
dimensions  of  Equation  2.2-5  and  for  plane  problems  the  terms  in  g 
corresponding  to  are  zero. 

The  element  stiffness  matrix  is  defined  as 

1  1 

K  ■  /bT  CB  -  /  /bTCB  I  det 

area 

where  C  is  defined  in‘  Equation  2.1-1  and  det  J 
Jacobian  of  the  coordinate  transformation.  The 
using  Gaussian  quadrature  as 


K  = 


where  W.  is  the  weight 

2.2.1  12-Node  Quadrilateral  Element 

The  element  is  shown  in  both  the  global  coordinate  and  its  mapped 
coordinate  systems  in  Figure  2. 2. 1-1.  The  edges  of  the  element  correspond 
to  values  of  s  or  t  of±  1,  and  the  midside  nodes  correspond  to  values  of  s 
or  t  of  -1/3. 

The  order  of  the  element  is  cubic  meaning  that  the  geometry  and 
displacement  components  can  vary  along  the  element  edge  as  a  third  order 
polynomial.  Strains  given  by  the  derivatives  of  the  shape  functions  may 
vary  quadratically  over  the  element.  The  specific  shape  functions  for  this 
element  are: 


IlBTCB|detJ| 

i  j 


l¥j 


V4j 


(2.2-7) 


function  corresponding  to  the  ith  quadrature  point. 


J  ds  dt  (2.2-6) 


is  the  determinate  of  the 
integration  is  performed 
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Global  Coordinate  System 


napped  Coordinate  System 

Figure  2.2. 1-1  12-Node  Quadrilateral  Element  and  Local  Element  Node  Numbers  . 
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=  32  0  “  t)(l  -  s)[-10  +  9  (s2  +  t2)] 


(2. 2. 1-1) 


N2  =  12  (1  '  t)(1  "  s2)0  ‘  3s) 

n3  *  §2  0  -  t)(l  -  s2)(l  +  3s) 

n4  5  1  (]  '  *)(1  +s)[-10  +  9  (s2  +  t2)] 

N5  =  32  (1  +  S>H  •  t2)0  -  3t) 

N6  3  33  (1  +  s)(1  *  t2^l  +3t) 

^7  =  32  (1  +  s)(l  *  t)[-10  +  9(s2  +  t2)] 
N8  =  33  (1  +  -  s2)(1  +  3s) 

Ng  *  ~  (1  +  t)(l  -  s2)(l  -  3s) 

N10  =  33  (1  +  -  s)[-10  +  9(s2  +  t2)] 

NH  =  33  ^  “  s)0  -  t2)(1  +  3t) 

N12  3  33  0  -  s)0  -  t2)(1  -  3t) 
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The  isoparametric  elements  can  be  "enriched"  to  account  for  the 
singularity  at  the  tip  of  a  crack  as  is  discussed  in  Section  2. 6.1.1.  The 
12-node  enriched  element  involves  the  addition  of  another  shape  function,  or 
creation  of  a  sub-paranetric  element.  The  functions  used  to  describe  the 
geometry  are  the  same  as  the  isoparametric  case,  but  the  displacement  field 
generated  by  the  crack  tip  singularity  is  added  to  the  displacements. 
Notation  for  this  will  be: 


(if  needed) 


(2. 2. 1-2) 


where  the  definition  of  these  fields  are  detailed  in  Section  2.6.1.  These 
modifications  to  the  shape  functions  are  also  reflected  in  the  element 
stiffness  matrix: 

B^B  +  B*  (2. 2. 1-3) 


K  =  f  (B  )TCB  da  (2. 2.1-4) 

area 


2.2.2  9-Node  Triangular  Element 

This  element  is  developed  as  an  extreme  case  of  the  12-node 
quadrilateral  elanent.  The  shape  functions  for  the  fourth  side  are  combined 
into  the  first  node.  Gaussian  quadrature  is  still  used  except  that  the 
accuracy  is  biased  by  the  skewing  of  the  integration  toward  the  first  corner 
of  the  element.  When  the  element  stresses  are  smoothed  to  the  nodes,  the 
sane  smoothing  is  applied  as  the  quadrilateral  except  that  the  stresses  for 
the  fourth  side  are  averaged  at  the  first  node. 

Since  we  use  the  collapsed  quadrilateral,  no  change  in  the  enrichment 
procedure  is  required. 

2.2.3  3-Node  Distributed  spring 

The  3-node  distributed  spring  was  developed  to  provide  the  same  order 
of  accuracy  in  elastic  boundary  conditions  as  is  provided  by  the  12-node 
quadrilateral  element.  The  stiffness  matrix  for  the  distributed  spring  is 
calculated  using  a  modified  line  integral  based  on  Gaussian  quadrature. 
Inputs  to  the  subroutine  are  the  axial  and  shear  stiffnesses  of  the  spring 


>r>N3dX3  1  N3WNM1  ac**>  iw  n-J-mnnuj-JM 
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Y  Global 


X  Global 


Figure  2. 2. 3-1  8-Node  Distributed  Spring 


Figure  2. 2.4-1  2-Node  Concentrated  Spring 
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defined  normal  to  the  first  side.  The  subroutine  converts  these  local 
values  into  global  stiffnesses.  Ttie  normal  at  the  node  point  makes  an  angle 
with  the  global  y-axis  of: 


a  =  Arctan 


(2. 2.3-1) 


,s'  refers  to  the  partial  derivative  with  respect  to  s. 


The  global  stiffness  components  at  a  node  become: 


KGx  =  C0S  a  *Lx 


2  a  K,  ..  +  Sin2  a  K, 


(2. 2. 3-2) 


KQy  =  Sin2  a  KLx  +  cos2  a  KLy 


(2. 2. 3-3) 


Kg  =  Sin  a  COS  a  (K[y  ■  Ku) 


(2. 2.3-4) 


where: 


=  local  axial  stiffness  at  that  node 


f<Ly  =  local  shear  stiffness  at  that  node 


Cnee  the  global  stiffnesses  are  known,  a  reduced  form  of  the 
integration  described  in  Equation  2.2-7  can  be  performed.  The  thickness  of 
the  element  in  the  t  coordinate  direction  does  not  enter  into  the 
formulation.  Therefore: 


K  =  J\dS 


(2. 2. 3-5) 


/*•>£  * 


=yKr  (s.). 

^  G . 


where  S  is  a  measure  of  position  along  the  spring  elanent  in  the  physical 
space  and  s  is  the  corresponding  position  in  the  mapped  space. 


2.2.4  2-Node  Concentrated  Spring 


The  2-node  concentrated  spring  allows  node-to-node  connections,  either 
axially  or  by  shear.  The  stiffness  matrix  is  determined  explicitly  from  the 
global  geometry. 
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In  local  coordinates  the  stiffness  matrix  is: 


where 

K  is  the  axial  stiffness 

Kx  is  the  shear  stiffness 
xy 

and  it  operates  in  the  nodal  displacement  vector  where  U  and  V  refer  to  the 
local  displacements  at  node  1  and  node  2  (see  Figure  2. 2. 4-1) 

The  stiffness  matrix  must  be  rotated  into  the  global  coordinate  system. 
This  is  done  explicitly  to  reduce  the  amount  of  formulation  time  required 
by  the  program: 


(2. 2. 4-2) 


where: 
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*  sin2  ~K 

+ 

cos2 

9K 

X 

xy 

*  the  angle  between  the  local  x  axis  and  global  x  axis 
2.2.5  Non-equilibration  of  Spring  Elements 

The  development  of  the  spring  elements  was  intended  to  supply  a 


/ti  Arthur  D.  Little,  Inc 


facility  for  modeling  of  both  axial  restraint  and  frictional  shear.  The 
axial  restraint  formulation  is  straightforward  and  widely  documented. 
However,  the  generality  of  the  formulation  vdiich  includes  shear  stiffness 
results  in  an  element  which  may  cause  non-equilibrating  moments  if  the 
element  has  a  finite  thickness  or  length.  Force  equilibrium  is,  however, 
maintained  unconditionally. 
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2.3  Nodal  Coordinate  Specification 


There  is  limited,  but  very  helpful,  mesh  generation  capabilities  in 
PAPST.  In  most  models,  only  the  corner  nodes  of  the  elements  are  required 
to  define  the  geometry.  Area  elements,  12  and  9  nodes,  can  be  specified 
with  just  4  or  3  nodes.  Intermediate  side  node  nunbers  and  coordinates  are 
automatically  generated  at  third  point  locations.  (Note:  This  implies 
straight  edges,  curved  edges  must  be  explicitly  specified.)  The  spring 
elements  also  have  generation  capabilities.  Distributed  springs  can  be 
specified  with  2  or  4  corner  nodes.  If  only  2  nodes  are  specified,  the 
program  will  generate  the  other  2  corner  nodes  at  the  same  global 
coordinates  and  ground  them.  Intermediate  nodes  are  generated  similar  to 
the  area  elements.  Concentrated  springs  can  be  specified  with  1  or  2 
nodes.  If  only  one-,  node  is  specified,  the  other  node  is  generated  at  the 
same  global  coordinates  and  grounded.  In  the  2  node  case,  both  nodes  are 
not  required  to  be  attached  to  the  area  elements. 

It  is  suggested  that  this  facility  be  used  vhenever  possible.  If  the 
user  decides  to  specify  the  intermediate  nodes,  the  third  point  spacing 
should  be  maintained.  Actual  coordinates  may  be  input  in  either  global 
rectangular  or  local  polar  coordinates.  The  progran  will  convert  all  local 
specification  to  polar  before  the  analysis  begins. 
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2.4  Nodal  Constraints 


Nodes  may  have  either  of  two  types  of  displacement  constraints,  the 
displacement  components  at  a  node  may  be  set  explicitly,  or  they  may  be  tied 
to  those  of  another  node.  Both  of  these  constraint  types  are  set  in  the 
local  coordinates  (user  specified)  of  the  node  or  nodes. 

2.4.1  Specified  Displacements 

The  value  of  either  or  both  dispacement  components  at  a  node  may  be 
explicitly  set  by  the  user.  The  following  information  refers  to 
implementation  of  displacement  specifications: 

1)  The  structure  must  be  constrained  against  rigid  body 
rotation  or  displacement. 

2)  Axisymmetric  problems  are  automatically  constrained  in  the 
x-direction  at  the  axis  of  symmetry.  Only  rigid  body  motion 
in  the  y-direction  must  be  constrained  by  the  user. 

3)  Constraints  may  be  specified  in  local  coordinates,  enabling 
the  user  to  model  inclined  supports. 

4)  Intermediate  side  nodes  will  be  constrained  automatically  in 
accordance  with  the  corner  node  constraints. 

5)  When  solving  symmetric  crack  problems,  a  synroetry  constraint  should 
be  applied  to  the  node  corresponding  to  the  crack  tip  as  well  as  to 
other  nodes  on  the  plane  of  symmetry.  Special  constraints  are 
required  if  the  user  if  modeling  tearing  (see  Section  2.6.4) 


2.4.2  Tied  Displacements 

Up  to  30  nodes  in  a  set  may  be  linked  together  to  have  common 
(but  unknown)  displacement  components .  Either  degree  of  freedom  at  a  node 
may  be  tied  to  that  at  another  node.  However,  each  degree  of  freedom  may  not 
be  tied  independently  to  separate  nodes.  Each  node  in  the  set  must  be  in  the 
same  local  coordinate  system. 
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2.5  Loadings 

;  Both  mechanical  and  thermal  loadings  are  available  to  the  user. 

1  Concentrated  forces  at  the  nodes  can  be  specified.  Distributed  tractions 

are  converted  to  consistent  nodal  forces  and  corrected  for  large 
5  displacement  effects.  Temperatures  at  each  of  the  nodes  can  be  included  in 

i  a  thermal  stress  analysis.  Pressures  applied  to  a  crack  face  include  the 

singularity  effect  in  the  consistent  nodal  forces. 

j  For  elastic  analysis  only,  all  loads  are  imposed  at  full  value. 

Thermal,  load  controlled,  and  displacement  controlled  analyses  are  all 
treated  separately  in  plastic  analyses.  Thermal  loadings,  mechanical  loads, 
•;  and  non-zero  displacements  cannot  be  handled  concurrently  in  the  plastic 

case.  The  program  will  handle  in  order  of  preference  thermal,  displacement, 
and  mechanical  loading,  and  ignore  additional  conflicting  information  if 
,  specified.  Load  factors  defined  in  the  non-linear  analysis  section  (2.6.2) 

:■  apply  to  the  temperatures,  the  forces  or  the  non-zero  displacements. 

2.5.1  Concentrated  Nodal  Loads 

l 

Concentrated  forces  can  be  applied  in  x  and  y  components  at  each  of  the 
unconstrained  nodal  degrees-of-freedom.  For  axisymmetric  problems, 

.  concentrated  nodal  loads  are  actually  line  loads  vhich  act  on  the  complete 

circunference  of  the  structure.  Such  loads  are  input  on  a  load  per  radian 
basis.  For  example,  if  a  load  of  3  units  per  circumferential  unit  of  length 
acts  at  a  distance  x  from  the  axis  of  symmetry,  the  total  force  F  is  3  (2~ 
x)  *  6rr  x.  The  required  input  load,  on  a  per  radian  basis,  is  F/2 r  or  3x. 

2.5.2  Distributed  Tractions 

Each  element  boundary  can  have  distributed  normal  and/or  shear 
tractions.  The  distribution  of  such  tractions  along  the  element  edge  can  be 
•  completely  general;  however,  most  problems  can  be  treated  with  a  constant  or 

linear  variation.  Normal  tractions  are  positive  when  directed  into  the 
element,  and  shearing  tractions  are  positive  vhen  acting  in  a  counter- 
|  clockwise  sense  along  the  element  edge.  The  program  automatically  computes 

the  equivalent  nodal  loads  for  applied  tractions.  Consistent  load 
generation  is  especially  important  in  large  deformation  calculations  since 
the  load  value  and  direction  will  change  as  the  structure  deforms. 

Pressure  loading  can  be  specified  along  an  element  edge  that  is  shared 
.  by  two  elements  within  the  structure.  The  sign  convention  in  this  case  will 

i  be  determined  by  the  direction  implied  in  the  nodal  input.  The  program  will 

"apply"  the  load  on  the  element  for  *hich  the  nodal  order  matches  the 
element  specification. 


2.5.3  Thermal  Loading 

For  a  thermal  loading,  temperature  information  must  be  provided  at 
least  at  each  of  the  corner  nodes.  If  the  temperatures  are  not  given  for 
some  (or  all)  pairs  of  intermediate  side  nodes,  PAPST  will  automatically  set 
those  values  from  a  linear  interpolation  of  the  values  at  the  corner  nodes. 
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In  certain  thermal  problems,  an  abrupt  change  in  temperature  may  exist 
at  an  interface  between  elements,  and  thermal  loading  is  desired  for 
elements  to  one  side  of  the  interface  but  not  for  those  on  the  other.  In 
this  case,  it  is  necessary  to  use  two  different  materials.  The  coefficient 
of  thermal  expansion  should  be  set  equal  to  zero  for  elements  for  which  no 
thermal  loading  should  be  calculated. 

PAPST  has  the  capability  to  do  thermal  loading  in  both  elastic  and 
plastic  analyses.  In  the  elastic  analysis,  the  thermal  loading  is  combined 
with  loads  from  pressures  and  imposed  displacements.  The  plastic  analysis, 
on  the  other  hand,  will  not  allow  any  other  loadings  to  be  applied  and  will 
override  them  if  specified  by  the  user.  All  of  the  non-linear  input 
specifications  remain  unchanged,  except  now  the  load  factors  are  applied  to 
the  temperatures. 

Thermal  loading  affects  two  sections  of  the  program,  calculation  of  the 
strain  components  and  body  forces. 

Thermally  induced  volumetric  expansion  must  be  subtracted  from  the 
strain  increment  in  order  to  calculate  the  mechanical  strain  increment: 


4  £  mech  .  4.  -  4«thermal 


(2. 5. 3-1) 


where  the  thermal  strains  are: 


A<thermal  .  alT 


(2. 5. 3-2) 
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Thermal  body  forces  are  calculated  vhen  the  element  stiffness  matrix  is 
formed.  These  nodal  forces  are  integrated  over  the  element  with  the  same 
Gaussian  quadrature: 

/•t  f.b./ii\jA  (2.5.3— 3) 

aFB  =  /Bt(u)C(u)a«  (u'dA 

area 

.  TX(btc  « f'b' (<iet 
i  j  I 


where  the  shape  functions  and  material  properties  are  at  the  current  values. 


'*D13/D33 

1  e  -  l-»23/033 

0 

_  "D43/D33 

for  plane 
strain 


The  D. .  terms  in  the  plane  strain  case  refer  to  the  terms  in  the  current 
material  property  matrix  before  they  are  modified  for  plane  strain  (see 
section  2.1.1). 
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2.6  OPTIONAL  ANALYSES 


The  analyses  contained  in  this  section  are  special  purpose 
developments  to  better  model  linear  elastic  and  elastic-plastic  fracture 
mechanics  problems.  Specifically  covered  are:  stress  intensity 
factor  calculation  (Reference  8) ,  J-integral  (Reference  1)  calculation,  and 
modeling  of  tearing  (Reference  2). 

2.6.1  Fracture  Mechanics 


The  concept  of  enriching  conventional  elements  with  the  asymptotic 
displacement  field  appropriate  for  fracture  mechanics  originated  with 
Benzley.  Extensive  application  has  been  made  for  both  Mode  I  and  Mode  II 
in  the  elastic  case.  Implementation  of  the  elastoplastic  case  is 
restricted  to  use  in  Mode  I  with  the  multilinear  material  model. 

The  approach  is  to  add  the  leading  terms  in  the  asymptotic  expansion 
to  polynomial  displacement  field.  This  adds  two  extra  unknowns  to  the  set 
of  nodal  displacements  of  which  one  is  zero  for  Mode  I  problems. 

The  asymptotic  crack  tip  singular  solution  is  described  in  the  next 
section.  The  implementation  of  the  combined  mode  enrichment  is 
subsequently  described  with  the  Mode  I  formulation  considered  as  a  special 
case. 

2. 6. 1.1  Asymptotic  Crack  Tip 

The  asymptotic  displacement  components  at  a  crack  tip  are  given  in  the 
standard  polar  coordinate  system  centered  at  the  crack  tip. 


3d 

u  =  VR/2tt  [5  -  3v  +  I2-  -  8b)cos  J  -  (1  +V  +  2  )cos|~]  (2. 6. 1.1-1) 

V  =  Ig-  yjR/ 2rr  [(7  -  V  +  _  8s)sine/2  -  (1  +  v  +  |^)sin  |^] 

For  Mode  II:  (KjjOnly) 

U  =  |g-  Vr/2T[(9  +  v  -  86)  sin  6/2  +  (1  +  v)  sin  39/2]  (2.6.1. 1-2) 

V  =  |_  yjR/z7  [-(3  -  5v  -  8e)  cos  e/2  -  (1  +  v)  cos  3e/2] 
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(2. 6. 1.1-3) 


1 


E  Ae-Acr 

slope  of  j-z  curve  (  — — — ) 
if  plastic  Kp 
0  if  elastic 


10  if  plane  stress 
(v  +  a/2)2/ (1  +  a)  otherwise 


(2.6. 1.1-4) 


with  the  angle  measured  from  the  crack  line  extension  path. 


2. 6. 1.2  Enriched  Elements 


The  enriched  quadrilateral  element  (Reference  9)  has  the  same  geometry 
as  the  standard  elements.  The  displacement  assumption  used  in  the  enriched 
element  is  essentially  the  usual  displacement  function  plus  the  leading 
terms  of  the  singular  displacement  expansion  in  Equation  2. 6.1. 1-1-2.  It 
takes  the  form 


2  2  3  7 

U  =  a-j  +  012s  +  +  a^S  +  a^St  +  aigt  +  a^S  +  agSt  (2. 6. 1.2-1) 

+  agSt  +  a-jQt  +  ojjSt  a-i2^s  +  ^i^i  (s,t)  +  Kjjg^(s,t) 


where  a.,  are  undetermined  coefficients 

Kj,  Kjj  are  the  stress  intensities  for  Mode  I  and  Mode  II  respectively 

f. (s,t)  g.  (s,t)  are  based  on  the  singularity  equations  (2. 6. 1.1-1, -2) 
and  1  1 

evaluated  in  terms  of  the  local  mapped  coordinates.  A  similiar  expression 
exists  for  the  v-component  of  displacement.  In  matrix  form,  Equation 
2. 6. 1.2-1  may  be  written  as 


U  (s,t)  =  P(s.t)  +  x,f,  (S.t)  +  Kng,  (s.t)  (2.6.1.2-21 
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3y  evaluating  Equation  2. 6. 1.2-2  at  each  of  the  nodes,  the  following 
matrix  equation  may  be  written 


U  =  O  +  *ifl  +  KII91 

(2. 6. 1.2-3) 


in  which  all  matrices  are  known  except  (J  and  a . 
Solving  for  a  then  gives 


a  =  C-1U  -  C”1  Kqf-j  -  C"1Kng1 
Substituting  for  a  then  gives 


(2. 6. 1.2-4) 


U(s,t)  =  PC~ 1 U  -  PC-1^^  -  PC'1Kng1  +  Kjf-j (s ,t)  +  Kj j 9 -j (s ,t) 


(2.6.1. 2-5) 


But  as  shown  in  Zienkiewicz,  (Reference  37,  page  106)  the  matrix  of 
displacement  interpolation  functions  i.e., 

1x12  12x12  1x12 

P  c  _1  =  N  =  [N1  N2  N3  .  .  .  .  N12]  (2. 6. 1.2-6) 


Consequently,  Equation  2. 6. 1.2-5  may  be  written 


u(s,t)  =  z  N.  u.  +  KT[f-,(s,t)  -  I  N.f,  ]  +  Kjj[g^(s,t)  -  -  N,.g^]  (2.6. 1.2-/) 

I  i  i  *  *  1  i  i 


where  the  second  level  subscripts  indicate  "evaluated  at  Node  i."  The 
analogous  expression  for  the  v-component  of  displacement  is 


v(s ,t)  =  Z  N.  V.  +  Kj[f2(s,t) 


(2. 6. 1.2-9) 


Equations  2. 6. 1.2-7  and  -8  are  the  displacement  approximations  used  in 
the  enriched  element.  It  should  be  noted  here,  however,  that  in  order  to 
provide  displacement  continuity  between  the  two  element  types,  the  second 
and  third  terms  on  the  right-hand  side  should  be  multiplied  by  functions  of 
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s  and  t  vdiich  cause  them  to  vanish  on  the  boundary  between  enriched 
elements  and  conventional  elements.  tote  also  that  the  functions  modifying 
the  stress  intensity  factors  multiplying  f^(s,t),  f^Cs^),  g,(s,t)  and 
g_(s,t)  lead  to  functions  of  s  and  t  in  the  26  x  25'<ielement  stiffness 
matrix  khich  cannot  be  numerically  integrated  with  the  same  precision  as 
the  functions  in  the  normal  stiffness  matrix.  Specifical Ly,  3  x  3  or  4  x  4 
Gaussian  quadrature  integration  is  adequate  for  the  conventional  element 
whereas  higher  quadrature  is  required  for  accurate  numerical  integration  of 
the  enriched  element  and  8x8  integration  is  used.  Higher  order 
integration  is  not  used  for  enriched  elements  in  the  nonlinear  program. 

The  significantly  snaller  enriched  element  used  in  the  plastic  analysis 
appears  to  obtain  acceptable  results  using  the  standard  4x4  quadrature. 

i 

2. 6. 1.3  Additional  Loading  From  Enrichment 

The  incorporation  of  the  singularity  displacement  field  has  an  effect 
on  the  loading  at  the  crack  tip.  3oth  the  distributed  pressure  and  the 
thermal  body  forces  must  be  modified  to  account  for  this  added  degree  of 
freedom. 

The  pressure  contribution  takes  the  form: 


r 

P  =  J*pUdr 

®  r  a, 

=  r/2  J  plids 

=  r/2  £  Wj  (2  Nip.)j  U. 

j 


V  £  •  o  «  ±  • 


3-D 


when  reduced  to  the  Gaussian  quadrature  format,  where: 


(2. 6. 1.3-2) 


4(1+  ct-  3 ) 

r  2t 


1 


where  S.  is  the  local  nodal  coordinate  of  node  j,  and  r  is  the  lencrth  of  the 
element^ edge  to  vhich  the  pressure  is  applied  and  the  other  variables  are 
consistent  with  the  notation  in  Section  2. 6. 1.1. 
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The  theraal  body  force  term  at  the  crack  tip  is  of  the  form: 


(2.5. 1.3-3) 

•-FB*  =  /  B*T(U)C(U)Ae  f,b-dA 
area 

where  the  notation  is  consistent  with  Section  2.5.3  and  2. 6. 1.2. 

2. 6. 1.4  Stress  Intensity  Factor  Recovery 

The  solution  of  the  set  of  governing  equations  for  the  global 
displacements  also  gives  the  normalized  stress  intensity  factors.  Recovery 
of  the  actual  stress  intensity  factors  becomes  trivial: 

Kj  =  Ex  Uj  (2. 5. 1.4-1) 


(2. 5. 1.4-2) 


Having  the  stress  intensity  factors  enable  us  to  calculate  the  strain 
energy  release  rate  at  the  crack  tip  (G  in  elasticity  is  comparable  to  J  in 
plasticity) .  The  following  equations  are  for  cracks  in  orthotropic 
materials  that  are  aligned  with  the  principal  material  axes.  Since  the 
program  will  only  accept  material  properties  in  global  coordinates,  this 
restricts  the  orientations  to  0°  and  90°.  The  equations  reduce  to 
isotropic  unconditionally  (Reference  38) . 
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2.6.2  J-Inteqral 


The  J-integral  (Reference  2)  is  a  measure  of  the  amplitude  of  the  near 
crack  tip  stress  and  strain  fields  in  an  elastic-plastic  material.  It  has 
been  proposed  as  a  fracture  initiation  criterion.  Under  conditions  of 
small  scale  yielding  for  which  linear  elastic  fracture  mechanics  is 
applicable,  J  reduces  to  the  energy  release  rate  G.  It  is  related  to  the 
Mode  I  stress  intensity  factor  by 


J  = 


i 

' 

< 

i 


Kj/E 

0-v2)  Kj/E 


for  plane  stress 
for  plane  strain 


(2. 6. 2-1) 


The  J-integral  for  elastic-plastic  crack  problems  is  given  in  terms  of 
the  strain  energy  density: 


(2.6. 2-2) 


W 


ij 


as 


J  = 


(2. 6. 2-3) 


where  the  crack  lies  along  the  local  x  axis  and  the  contour  is  traversed 
counterclockwise  from  one  crack  face  to  the  other.  Figure  2. 6. 2-1  shows 
the  general  contour.  The  traction  vector  has  components  T.  and  the 
corresponding  displacement  components  are  U^.  The  differential  length 
along  the  arc  is  indicated  by  ds. 

In  the  case  of  finite  elements,  a  path  must  be  traversed  along  the 
element  boundaries.  Typical  path  contours  are  shown  in  Figure  2. 6. 2-1. 

The  integral  is  performed  using  the  averaged  nodal  stresses  and  strains  and 
a  line  Gaussian  quadrature  integration.  The  use  of  averaged  values  avoids 
any  element  or  path  dependence. 


The  J-integral  is  calculated  as  follows: 


.  35 
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(2.S.2-4) 


J  = 


S  J[s1n  5 2laxy  s+ Vy  ■  w) 


edges  s 


+  cos  v>x)]  dx 


2  xy  x  y 

+  [cose2(W  -  cx£x  -  Jxy  V,x  ) 

-sins2(cxU’y  +  cxyey)-1  dy 


where:  (for  power  hardening  materials) 


i  (ii+jiL*  +-iT,(  (^-)n+1  -  n°Ll 

E  L  3  e  6  PP  n  +  1  Gyd  yd  J 


(2. 6. 2-5) 


W  =  i 


(for  multilinear  materials) 


m-1 

LL ±v)a  2  +  (1_-  2v^2  +  y  !i  (  2 
3  ce  +  6  "on  +  Zf 


'1+1 


iy  .  -n  /  2  2x 
°i}  +  2~  (oe  *  °m} 


and: 


cr  „  ~  g y v  +  c  + ' a 

PP  xx  yy  22 


v  =  v 


xy 


(2. 6. 2-5) 


(2. 6. 2-7) 


E  =  E 


An  additional  contribution  must  be  included  for  the  axisymmetric 
j-integral  (Reference  10) : 


(2. 6. 2-8) 

Jarea  =  _AC0S92<V2  '  Vx  ’  °xy  V-x>/z 

A 

+  sin«n(a  U,,.  -  c  c  )/z]  dA 
?  X  y  xy  y"  J 

This  area  integral  is  added  to  the  contour  integral  for  all  elements 
that  fall  within  the  path  contour,  i.e.. 


J  *  J  . .  +  J 

path  area 


(2. 6. 2-9) 
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The  energy  release  rate  defined  in  Section  2. 6. 1.4  can  use  the  plastic 
singularity  to  approximate  the  J  value  for  a  multilinear  material  model. 

If  is  substituted  for  K_  in  equation  2. 6. 1.4-3,  G  becomes  J  with  the 
following  modifications:  1  ° 


For  plane  stress: 


Jp  =  0  +  a)  Gj 


(2.6.2-10) 


For  plane  strain: 

[1  +  y^/2  ~  3vy/2  +  am(  1  -  3y/4  +  y*V2)1 

J  =  n  2.  G  (2.5.2-11) 

P  U  “  v  )  1 

where: 


(2.6.2-12) 


Y  =  (v  +  an/2)/(l 


and  c^is  the  slope  of  the  last  ligament  of  the  plastic  stress-strain 
relationship  defined  by  equation  2. 1.4-1 


2.6.3  Tearing  Modulus 


The  applied  Tearing  Modulus 


14 


is  defined  by 


TAPPL 


E  dJ  _  E  M 

ZZ  da  ~  a  2  Aa 

o  0 


(2. 6. 3-1) 
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E  is  Young's  modulus,  o,  is  the  flow  stress  (average  of  yield  and  ultimate 
stresses) ,  and  dJ/da  represents  the  rate  of  change  of  the  J-integral  with 
respect  to  crack  length.  Ductile  crack  growth  stability  is  assessed  using 
the  tearing  modulus  as  follows:  Fracture  is  stable  upon  reaching  if 
<  TvaT.,  or  unstable  if  T  L  >  T  w»mT ,  where  T 
oauius  of  the  material  taxen  from 


‘APPL 

tearing  mod 
of  J  versus  Aa  in  standard  JIC  tests, 


’t  c 
is  th§^ 


.  '  ’"“‘j  *  .MATL  t 

est  records  (resistance  curves) 


Finite  element  analysis  including  crack  growth  (Reference  11)  can  be 
performed  to  calculate  T.  In  PAPST,  crack  extension  over  one  element  can 
be  analyzed.  This  Js  done  by  inserting  a  very  stiff  distributed  spring 
along  the  path  defined  for  crack  growth.  The  user  is  allowed  to  define  the 
JIC  value  at  'fchich.  the  crack  should  grow.  Cnee  this  value  is  reached,  the 
spring  is  progressively  softened  holding  all  other  boundary  conditions 
fixed.  When  the  program  converges,  a  new  value  for  J  is  computed.  The 
Tearing  Modulus  can  then  be  calculated  from  Equation  2. 6.3-1. 
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2.7  OPTIONAL  FEATURES 


Several  optional  features  are  available  in  the  program  to  aid  the  user 
in  describing  the  model  and  interpreting  the  results.  This  section 
describes  the  user  selectable  features. 

2.7.1  Plotting 

The  plotting  package  has  the  following  capabilities: 

a)  Drawing  deformed  shape 

b)  Stress  contour  plots 

c)  Thermal  contour  plots 

d)  Enlargement  of  Selected  Regions 

This  package  has  been  developed  for  use  on  a  CALCCMP  plotter. 

2.7.2  Wavefront  Reordering 

This  feature  is  still  in  development.  The  user  can,  however,  utilize 
the  manual  reordering  feature.  This  allows  the  frontal  solution  to  be 
performed  in  a  different  order  from  the  element  input  ordering.  See 
Appendix  A,  Topic  6,  Subject  C  for  details. 
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3.0  VERIFICATION 


Verification  of  a  computer  program  is  an  on-going  process.  Each  new 
feature  or  modification  to  the  program  must  be  shown  to  work  itself  and  be 
shown  bo  not  affect  the  performance  of  other  existing  features  of  the 
program.  Considering  the  amount  of  modification  that  has  been  made  to 
PAPST  during  this  contract,  it  was  appropriate  to  collect  and  review  the 
previous  verificatioon  cases.  Fran  this  collection,  selected  test  cases 
were  used  to  baseline  this  version  against  previous  versions  to  see  what 
differences,  if  any,  existed,  and  if  they  did,  were  they  a  result  of 
improvements  and  modifications.  New  features,  such  as  springs  and  tearing, 
required  verification  to  theoretical  or  experimental  results. 

This  verification  effort  was  not  exhaustive.  There  are  many  features 
and  combinations  of  features  that  have  not  been  explicitly  verified. 
Experience  with  the  program,  both  here  and  at  the  Navy,  leads  us  to  believe 
that  there  are  no  obvious  problems  or  errors. 

3.1  General  Tbpics 

3.1.1  Elastic  Analysis 

This  version  of  PAPST  was  developed  from  the  most  recent  version  of 
APES.  The  broad  base  of  experience  with  the  use  of  APES,  which  is  now 
fully  incorporated,  implies  that  this  section  of  the  program  is  partially 
verified.  Also,  efforts  required  in  later  portions  of  this  section  require 
that  the  elastic  solution  be  correct.  Therefore,  no  specific  effort  was 
undertaken  to  verify  this  topic. 

3.1.2  Elastic  Singularities 

In  Section  3.5  one  of  the  comparisons  to  verify  the  j-integral 
involved  the  use  of  the  elastic  singular  solution.  The  results  from  the 
elastic  singular  solution  were  found  to  be  consistent  with  previously 
reported  values. 

3.1.3  Multiple  Runs 

Several  of  the  verification  test  cases  were  run  sequentially  through 
the  program.  Included  were  multiple  elastic  and  multiple  elastic/plastic 
cases  as  well  as  combinations  of  elastic  and  elastic/plastic  cases.  All  of 
the  problems  ran  identical  to  cases  run  individually. 

3.1.4  Inclined  Supports 

Rotated  test  cases,  similar  to  the  J „  test  cases  used  in  Section  3.5, 
were  generated  and  run  to  verify  that  the Solution  was  insensitive  to  the 
model  orientation.  Table  3.1-1  summarizes  the  results.  Overall  comparison 
is  excellent.  The  only  noticeable  discrepancy  arises  fran  the  plastic 
singularity,  and  this  vras  expected  because  it  is  not  generally  accurate. 
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TABLE  3.1-1 


COMPARISON  OF  INCLINED  MODEL 

TO  UNROTATED 

MODEL 

DESCRIPTION 

0=0° 

0=45° 

A% 

Model  #1  - 

6EL,  Plastic,  Plane  Strain 

Incr.  ?1 

PFRAC  (1) 

1.312940 

1.312767 

Highest  Loaded  Node 

62038.339 

62040.306 

- 

o1 

J2 

85.010 

85.003 

- 

65.263 

65.256 

.011 

Incr.  #2 

Highest  Loaded  Node 

110050.114 

111451.825 

1.27 

J2  , 

803.70 

776.55 

3.378 

696.78 

690.79 

0.86 

Incr.  #3 

Highest  Loaded  Node 

139335.383 

140941.39 

1.15 

2079.0 

2048.7 

1 .46 

1690.06 

1682.8 

0.43 

Model  #2  - 

6EL,  Plastic,  Plain  Strain, 
Enriched 

Incr.  #1 

PFRAC  (1) 

1.244933 

1.244870 

K_ 

5103.45 

5703.87 

- 

o? 

H?ghest  Loaded  Node 

77.103 

77.108 

- 

62.876 

62.878 

- 

85189.317 

85210.081 

.02 

Incr.  #2 

kd 

28066.948 

30595.377 

9.01 

Highest  Loaded  Node 

130250.355 

133047.625 

2.15 

J2 

803.75 

775.62 

3.5 

689.61 

683.47 

0.89 

Incr.  #3 

kd 

57500.24 

60063.13 

4.46 

Highest  Loaded  Node 

174047.27 

176820.054 

1.59 

o’ 

2066.7 

2033.5 

1.61 

1743.8 

174.06 

0.81 

*  <.01% 
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3.2  Concentrated  Springs 

The  concentrated  spring  verification  addresses  each  of  the  possible 
uses  of  the  spring.  Features  of  the  spring  are: 

1)  Axial  and  shear  stiffness 

2)  1-node  -  grounded  (zero  length) 

3)  2-node  -  between  element 

4)  Rotated  or  inclined  specification 

5)  Combined  stiffness  action 

6)  Potential  violation  of  moment  equilibrium 

« 

Four  test  cases  were  developed  to  cover  all  of  these  topics.  Three  of 
the  cases  are  statically  determinate,  permitting  easy  checking.  The  fourth 
case  demonstrates  the  problem  with  moment  equilibrium. 

Test  Case  #1,  shown  in  Figure  3.2-1,  is  one  element  that  is  pinned  at 
one  end  and  restrained  by  a  spring  at  the  other.  An  overturning  force  is 
applied.  The  spring  was  defined  by  one  node,  resulting  in  the  grounding  of 
the  other  end  and  zero  length.  Cnly  the  shear  stiffness  was  input.  The 
local  coordinate  system  defaults  to  the  global  system  and  thereby  defines  a 
spring  with  vertical  stiffness  only. 

Test  Case  #2,  shown  in  Figure  3.2-1,  places  an  axial  spring  between  two 
elements.  Both  elements  are  restrained  from  vertical  motion.  The  left 
element  is  pinned  at  the  far  side  and  the  right  element  has  a  uniform 
traction  applied.  The  force  from  the  traction  must  be  transferred  through 
the  spring  to  the  far  support.  Stabilizing  moments  at  the  base  of  each 
element  were  checked  for  equilibrium. 

Test  Case  #3,  shown  in  Figure  3.2-2,  has  an  axial  spring  that  is 
defined  between  two  pinned  supports.  One  of  the  supports  undergoes  a  known 
displacement.  Resultant  spring  forces  are  axial. 

Test  Case  #4,  shown  in  Figure  3.2-2,  has  the  same  model  geometry  as 
Test  Case  #3.  The  difference  is  that  both  the  axial  and  shear  stiffnesses 
were  set  to  the  same  value.  This  results  in  an  equal  restraint  against 
motion  in  all  directions.  The  force  at  the  displaced  end  of  the  spring  is 
the  spring  stiffness  times  the  displacement.  No  other  forces  are  developed. 
The  other  end  of  the  spring  has  a  support  reaction  that  is  equal  and 
opposite  to  the  force  developed  at  the  displacement.  Therefore,  moment 
equilibrium  is  not  maintained  over  the  finite  length  of  the  spring. 

The  computer  runs  that  were  made  agree  completely  with  the  analyses 
shown  in  Figures  3.2-1  and  3.2-2. 
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s  * 


U  =  2P/K  , 
xy 


K  =0 
x 


Example  #1 


P=2pa 


Figure  3.2-1  Concentrated  Springs  -  Examples  #1  and  #2 
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t-  =  K  u  cos  e 
ax  « 

Fx  =  Kx  u  cos ^9 

Fy  =  Kx  u  sin  9  cos  9 


Example  #3  Axial  and  Rotated 


K*  -  V  ■ K 


M 


F,  K  \ 


V  '54 


Fx  =  K  u 

Fz  =  K  u  COS  9 

Fy  =  K  u  sin  e 


Example  #4  Combined  and  Rotated 
Figure  3.2-2  Concentrated  Springs  -  Examples  #3  and  #4 
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3.3  Distributed  Springs 


The  distributed  spring  verification  is  similar  to  that  for  the 
concentrated  spring.  Most  of  the  possible  uses  and  features  were  tested. 
However,  the  added  complexity  of  the  distributed  spring  makes  an  identical 
level  of  study  impractical.  Only  the  most  common  features  are  addressed, 
and  we  rely  on  the  generality  of  the  derivation  to  cover  other,  less  common 
cases . 


Features  of  the  springs  are: 

1)  Axial  and  shear  stiffness 

2)  2-node  -  grounded  (intermediate  nodes  generated) 

3)  4-node  -  between  element  (intermediate  nodes  generated) 

4)  Rotated  or  inclined  specification 

5)  Combined  stiffness  action 

6)  Potential  violation  of  moment  equilibrium 

Three  test  cases  were  developed  to  cover  these  topics.  The  one  case 
that  demonstrates  the  potential  violation  of  moment  equilibrium  is 
statically  determinate  if  this  is  kept  in  mind.  The  other  two  can  be 
quickly  analyzed  by  hand. 

Test  Case  #1,  shown  in  Figure  3.3-1,  is  a  row  of  elements  resting  on  an 
elastic  foundation  provided  by  the  springs.  Uniform  suction  is  applied  to 
the  top  surface.  The  springs  are  defined  by  2  nodes  and,  therefore,  have 
zero  length  and  are  grounded.  Axial  stiffness  in  this  case  corresponds  to  a 
vertical  restraint.  The  distribution  of  pressure  on  the  top  surface  is 
reproduced  in  the  reactions  provided  by  the  springs. 

Test  Case  #2,  shown  in  Figure  3.3-2,  is  the  same  general  configuration 
as  Test  Case  #1.  Rather  than  suction,  uniform  shear  is  applied  to  the  top 
surface.  Vertical  restraints  are  defined  along  the  bottom  edge  of  the 
elements.  The  springs  have  a  finite  length  and  only  a  shear  stiffness. 
Deformation  of  the  elements  causes  some  small  redistribution  of  shear  force 
to  the  springs.  The  moment  developed  by  the  vertical  restraint 
counterbalances  the  moment  caused  by  the  shear  force.  However,  this  moment 
has  a  magnitude  that  implies  that  the  springs  are  supplying  the  reaction  at 
the  point  of  connection  to  the  element  and  not  at  the  ground  of  the  springs. 
This  is  consistent  with  the  warning  on  the  violation  of  moment  equilibrium. 
The  formulation  of  the  spring  does  not  account  for  length,  and  therefore, 
shear  springs  of  finite  length  will  not  balance  properly. 

Test  Case  #3,  shown  in  Figure  3.3-3,  demonstrates  the  use  of  the  spring 
between  elements.  Two  elements  are  connected  by  a  spring  on  an  inclined 
surface.  Suction  is  applied  to  the  unsupported  element  and  must  be 
transferred  through  the  spring  to  the  other  element  to  reach  a  support.  The 
moment  developed  at  the  support  balances  the  moment  from  the  suction.  Equal 
axial  and  shear  stiffnesses  provide  a  uniform  stable  connection  between  the 
two  elements.  A  biased  connection  could  be  specified  by  the  use  of 
different  spring  stiffnesses  for  axial  and  shear  components. 

Computer  runs  made  for  these  cases  agree  completely  with  the  analyses 
shown  in  Figures  3.3-1  through  3.3-3. 
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7.5 


Spring  Reactions 


Figure  3.3-2  Distributed  Springs  -  Example  #2 
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3.3-3  Distributed  Springs  -  Example  #3 
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3.4  Thick  Walled  Cylinder 

The  plane  strain  idealization  of  an  internally  pressurized  thick  walled 
cylinder  (an  axisymmetric  problem)  is  shown  in  Figure  3.4-1.  This  problem 
was  employed  to  verify  the  previous  version  of  PAPST.  The  method  of  load 
incrementation,  convergence  characteristics,  displacements  at  the  inside  and 
outside  surfaces  are  all  indicated  in  Thble  3.4-1  and  comparison  is  made 
with  the  results  of  Reference  12.  We  have  added  the  results  from  the 
present  PAPST  program  to  Table  3.4-1. 

In  order  to  also  correlate  these  results  to  analytic  example,  we 
limited  the  analysis  to  the  small  scale  displacement  option.  The  outer 
displacement  of  an  elastic/perfectly  plastic  cylinder  can  be  approximated 
for  a  Von  Mises  yield  criteria  from  work  done  by  Hill  (Reference  13) .  The 
elastic/plastic  boundary  can  be  found  by  the  equation. 


P/yd  =  ln(c/a)  +  V2(l  -  c2/b2) 


where : 

a  -  internal  radius 
b  =  external  radius 

c  =  elastic/plastic  interface  radius  a  c  b 
p  =  internal  pressure 

Jyd  =  stress  x  2/  /3~  to  account  for  Von  Mises  criterion 

rather  than  Tresca. 


The  displacement  field  in  the  elastic  region  then  becomes: 


U  =  a  yd  T?  (  0 -2v)r  +  b2/r  ) 

J  b 

where: 

E  =  Young's  modulus 
v  =  Poisson's  ratio 
r  =  radius  to  be  evaluated 


(3.4-2) 


Results  for  this  equation  are  compared  to  the  computed  results  on  Table 
3.4-1. 
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(12-Ncde  Eiements) 


2  345  678  9  10 

e  (Young's  Modulus)  =  6.895  (10)4  MPa 
v  (Poisson's  Ratio)  =  0.33 

oyci  (Yield  Stress)  «  68.95  MPa  Elastic  —  Perfectly 
'  "  Plastic  a  -  e  Curve 


Approximate  Locations  of  Integration  Points  at  which  Stress-Strain 
History  is  Monitored  are  Shown  by  Dots 


Figure  3.4-1  Idealization  of  Thick  Wall  Cylinder  for  Plastic  Analysis 


TABLE  3.4-1 


RESULTS  FOR  3  ELEMENT  IDEALIZATION  OF 
_ THICK  WALLED  CYLINDER 


Internal 


Increment 

Number 

Stress 

Small  Scale  Disc! . 

Larqe  Scale  Disp. 

Yield 

Stress 

Theory 
(Ref.  13) 

Initial 
Stress 
(Ref  .17) 

New  Papst 

Old 

PAPST 
(Ref.  18) 

1 

0.453 

NA 

1 

1 

.261 

.261 

.1616 

.161 

.161 

Wm&DSm 

2 

0.5 

3 

2 

2 

2 

.294 

.294 

.294 

.294 

.1811 

.181 

.181 

.181 

.181 

3 

0.55 

3 

2 

2 

2 

.332 

.337 

.337 

.338 

.2053 

.203 

.206 

.206 

.206 

4 

0.6 

4 

2 

2 

2 

.391 

.387 

.387 

.388 

.2348 

.236 

.234 

.235 

.235 

5 

0.65 

4 

2 

2 

2 

.455 

.456 

.456 

.457 

.2723 

.272 

.272 

.272 

.273 

6 

0.7 

5 

3 

3 

3 

.546 

.546 

.547 

.549 

.3215 

.321 

.321 

.321 

.323 

7 

0.75 

6 

3 

3 

3 

.688 

.689 

.691 

.696 

.3951 

.395 

.396 

.397 

3 

0.77 

7 

3 

3 

3 

.770 

.775 

.781 

.787 

.4391 

.438 

.440 

.444 

.447 

9 

0.79 

8 

3 

3 

3 

.898 

.923 

.936 

.950 

.5094 

.504 

.516 

.523 

.531 

(1)  Number  of  iterations  to  convergence  2 

(2)  Radial  displacement  at  inside  surface,  in  x  10  2 

(3)  Radial  displacement  at  outside  surface,  in  x  10^ 
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A  limited  amount  of  verification  of  the  option  has  been  performed  to 
assess  the  validity  of  the  J  integral  calculation.  The  previous  revision  of 
PAPST  had  extensive  testing  and  was  subjected  to  correlation  with  alternate 
solutions  (Reference  14).  While  the  calculation  procedure  has  been 
completely  rewritten,  the  basic  equations  and  approach  remain  unchanged. 

A  single  model  was  chosen  to  compare  against  previous  results.  The 
geometry  of  the  model  is  shown  in  Figure  3.5-1.  It  is  approximately  the 
same  as  the  geometry  used  for  Reference  14. 

I 

Three  separate  cases  were  analyzed:  elastic  enriched,  plastic,  and 
plastic  enriched.-'  The  cases  were  run  for  both  plane  strain  and  plane 
stress.  Results  for  these  runs  are  summarized  in  Thble  3.5-1.  Comparable 
results  are  also  shown  for  calculations  done  with  the  previous  revision 
of  Papst  (Reference  14) . 
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Oracle  Tip 


E  *  30  x  10^  psi 
v  =  .  3 


CTyd  =  psi 

a  =  10 
n  =  1 
s  =  0.0 


Plane  Stress 
and 

Plane  Strain 


Figure  3.5-1 
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3.6  Tearing  Modulus 


Three  separate  models  of  a  compact  tension  specimen  were  generated  to 
verify  the  tearing  calculation.  In  all  cases,  the  material  was  bilinear 
with  a  yield  stress  of  138  KSI  and  a  stress  level  of  152  KSI  at  10% 
elongation.  Young's  modulus  was  29,000  KSI,  Poisson's  ratio  was  .3,  and  a 
value  of  890  in-lb/in  was  'used  for  J_c.  These  properties  roughly 
correspond  to  HY-130. 

Each  of  the  models  contains  a  different  level  of  sophistication. 

Model  #1,  shown  in  Figure  3.5-1,  was  6  elements  with  the  crack  extension 
of  0.1  inches.  Model  #2,  shown  in  Figure  3.5-2,  was  18  elements  with  a 
crack  extension  of  0.05  inches.  Model  #3,  shown  in  Figure  3.6-3,  was  28 
elements  with  a  crack  extension  of  0.02  inches.  The  Model  #2  improvement 
over  Model  #1  was  to  improve  the  grid  in  the  vicinity  of  the  crack  tip. 
Model  #3  enhanced  this  grid  even  more  by  providing  a  refined  grid  at  the 
end  of  the  crack  growth  region.  It  also  refined  the  grid  in  the  region  of 
compression  yielding  near  the  back  surface. 

Each  of  the  models  was  loaded  in  four  increments.  Displacement 
control  was  imposed  on  a  spring  in  series  with  the  specimen.  Different 
spring  stiffnesses  were  tried  on  Models  #1  and  #2  to  approximate  changes 
in  machine  compliances. 

For  each  model  and  spring  combination,  two  runs  were  made.  The  model 
was  first  run  with  the  initial  crack  and  then  grown  to  the  final  length. 
Next  the  model  was  run  with  the  full  length  crack.  This  was  done  to 
compare  the  stress  state  of  the  growing  crack  to  the  full  length  crack  and 
to  check  for  any  numerical  anomalies  that  may  have  occurred.  None  were 
seen. 


The  J  value  that  was  used  for  comparison  was  the  path  that  fell 
entirely  within  the  elastic  region  of  the  specimen.  This  value  compares 
the  best  with  the  J  integral  estimate  based  on  the  applied  load-load  line 
displacement  curve.  At  this  time,  it  is  felt  that  the  deviation  of  the 
J-integral  path  values  (previously  reported)  is  a  result  of  the  anoothing 
technique  that  is  used  to  determine  the  nodal  stresses  and  strains.  While 
the  anoothing  technique  works  acceptably  for  elastic  strain  fields,  the 
complexity  that  arises  from  the  elastic-plastic  strain  fields  is  not  well 
represented  by  the  smoothing  functions. 

Figures  3.6-4  through  3.6-8  show  the  J  vs.  load  line  displacement 
curves  for  the  five  cases  reported.  It  can  be  seen  that  the  higher 
refinement  in  Model  #3  results  in  closer  comparison  between  the  growing 
and  the  grown  crack  predictions.  Figure  3.6-9  shows  the  plastic  zone 
comparison  between  them.  It  should  be  noted  that  the  strain  level  near 
the  crack  tip  was  about  30%  in  both  cases.  This  is  far  in  excess  of  the 
assumed  ultimate  of  10%. 

Figure  3.6-10  shows  the  tearing  modulus  values  that  were  calculated 
from  the  five  cases.  The  experimental  curves  that  are  used  for  comparison 
are  from  a  paper  by  Joyce  and  Vassilaros  (Reference  15) .  Several 
important  differences  should  be  noted: 
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1)  The  material  has  an  ultimate  stress  of  152  KSI  at  20% 
elongation.  The  finite  element  models  used  a  slightly  higher 
hardening  material . 

2 

2)  The  material  has  a  reported  J,c  of  873  in-lb/in  while  the  "key 
curve"  used  to  generate  the  T^p  curve  appears  to  have  a  J  of 
about  800  in-lb/in  .  The  models  all  grew  cracks  near  893 
in-lb/in  . 

3)  The  tearing  curves  are  for  continuously  growing  cracks.  The 
models,  by  necessity,  can  only  grow  the  cracks  incrementally. 

i 

4)  The  model  is  idealized  as  plane  stress.  The  specimens  are 
finite  thickness. 


£&N3dX3  INIWNMIAnt  <  W 
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Plane  Stress 
E  =  29  x  106  psi 
v  =  .  3 

ay(j  =  138  ksi 
a  =  152  ksi  at  e  =  10% 

3  *  l.G 

Figure  3.6-1  Model  #1-6  Element 
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Figure  3.6 
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Figure  3.6-8  Model  #3,  J  vs.  Load  Line  Displacement,  Km  =  60  K/i 
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Tearing  Modulus 


3 . 7  The  final  Load  i  ng 


For  verification  of  this  coding,  three  separate  models  were  each  run 
for  axisymmetry,  plane  strain,  and  plane  stress: 

1)  Unconstrained  with  a  uniform  temperature  field. 

2)  Constrained  top  and  bottom  with  a  linearly  varying  temperature 
field. 

3)  Constrained  top  and  bottom  with  a  uniform  temperature  field  and 
described  with  bilinear  elastic/plastic  material  properties. 

4 

Both  the  unconstrained  and  constrained  models  are  shown  in  Figure 
3.7-1.  In  all  cases.,  the  computer  model  matched  the  analytical  equations. 

3.7.1  Axisymnetry 


Test  case  one,  using  the  unconstrained  model  and  uniform  temperature, 
produced  uniform  strains  and  no  stresses  or  reactions.  This  is  the  trivial 
theoretical  case,  but  shows  that  no  residual  or  unbalanced  forces  are 
incurred  by  the  thermal  body  forces. 

Test  case  two  used  the  constrained  model  with  a  lineaerly  increasing 
temperature  field.  This  case  is  described  generally  in  Reference  16.  The 
stresses  in  the  body  are: 

(3.7-1) 


where  a  =  outer  radius 

r  *  a  radius  within  the  body 

T  *  temperature  field  expressed  as  a  function  of  r 
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Figure  3.7-1  Thermal  Models 
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"frie  linear  field  that  was  used  can  be  expressed  as: 


T  = 


T  r 
o  a 


(3.7-2) 


where  Tq  is  the  tsnperature  at  r  *  a 


Incorporating  Equation  3.7-2  into  3.7-1  and  integrating,  we  get  the 
following  stresses : 


aET 


a  = 


r  TO-v) 
aET. 


a.  = 


9  3(l-v) 


n  -  f] 


['•ri 


(3.7-3) 


.  _  aETo  r0  3r-, 


We  can  evaluate  the  term  in  brackets  at  the  center  and  at  the  outer 
radius: 


r=0  r=a 


o 


r 


a 


9 


1 

1 


0 

-1 


a 


Z 


2v  2v-3 


Test  case  three  also  used  the  constrained  model  but  with  a  uniform 
tsnperature  field.  The  elastic  stresses  in  this  case  become: 


°9  =  ° 
oz  =  -aET 


The  plastic  response  was  described  using  a  bilinear  material  model, 
(n  *  1).  The  stress  beyond  yield  in  this  case  becomes: 


*&N3<ixa 
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where: 

am  =  slope  of  the  stress/strain  curve  beyond  yield 

~i  =  tempereature  at  yield 

AT  =  temperature  increment  beyond  yield 

This  behavior  was  verified  in  the  computer  model. 
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3.7.2  Plane  Stress 


Plane  stress  is  quite  similar  to  axisymmetry  in  terms  of  the  stress 
fields  for  this  analysis.  Test  case  one,  unconstrained  with  uniform 
temperature,  produces  uniform  expansion  with  zero  stress.  Test  case  three 
with  the  same  material  properties  results  in  uniaxial  stress  of  the  form: 


(3.7-6) 

ay  =  "aETl  ‘  “ThT-  &T 
m 

Test  case  two  shows  different  behavior.  We  can  treat  the  case  of  a 
constrained  element  with  a  linearly  increasing  temperature  field  as 
unconstrained  expansion  with  a  linearly  increasing  imposed  displacement. 
This  displacement  can  be  described  as: 


(3.7-7) 


v(x)  =  -aaT(x)  =  -ctTQX 


The  total  strains  therefore  are: 


ey  =  v(x)/a  +  aT(x)  =  0 

e  =  -v(v(x)/a)  +  aT(x)  =  ( T+vJaTgX 
X 

e  =  -v(v(x)/a)  +  aT(x)  «  (l+v)aTQX 


(3.7-8) 


The  stress  field  becomes  therefore: 


a 


x 


=  0 


oy  =  -aET(x)  »  -  ctETQX 


(3.7-9) 


m 
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3.7.3  Plane  Strain 


Plane  strain  does  not  have  a  trivial  test  case  as  do  axisymmetry  and 
plane  stress.  The  Z  direction  is  always  constrained  and  stresses  are 
produced  in  that  direction  for  any  thermal  loading. 

Test  case  one  produces  the  following  total  strains  and  mechanical 
stresses : 


e  =  (l+v)aT  (3.7-13) 

A  U 

“y  =  (l+v)aT0 


Test  case  tv*5  can  be  analyzed  similarly  to  plane  stress,  except  that 
both  y  and  Z  directions  are  constrained.  The  mechanical  displacements 
required  to  offset  the  thermal  expansion  can  be  described  as: 


v(x)  =  -aaT(x) 
w(x)  =  -atT(x) 


(3.7-11) 


where  a  and  t  are  the  in-plane  and  transverse  dimensions  of  the  plate, 
respectively. 

The  only  non-zero  strain  is  in  the  x  direction. 


ex  =  frvaT^x)  +  cJ(x) 


(3.7-12) 


The  resulting  stress  state  is: 


ay  =  -aET(x)/(l-v) 
az  =  -oET(x)/(l-v) 


* 


(3.7-13) 
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When  these  are  evaluated  for  the  linear  temperature  distribution  in 
Equation  3.7-2,  they  check  with  the  results  of  the  computer  model. 

Test  case  three  results  in  the  same  strain  field  as  described  in 
Equation  3.7-12  except  that  the  temperature  field  is  constant  rather  than 
linear.  The  stresses  become  in  this  case: 


[ET.  + 


ok) 


aT] 


(3.7-14) 


where  the  notation  is  similar  to  Equation  3.7-5. 
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APPENDIX  A  -  CHANGES  IN  THE  INPUT  FILE  FOR  PAPST 


■  . 


1)  Crack  Growth  - 

Group  VIII  -  Nonlinear  Analysis  Data 
C.  J-Integral  Related  Data 

Card  3  -  Needed  only  if  NCNR  of  card  1  is  now  zero. 

The  change  in  the  crack  growth  mechanisn  required  a  modification  to  the  user 
input.  Only  one  increment  of  crack  growth  is  allowed.  The  procedure  is 
described  in  Section  2.6.3. 


CARD  3  -  NCN1,  NCN2*,  NON,  JIC,  SIGMA3 


FORMAT  (315/  2F10.2) 


(1-5)  NCN1  - 
(5-10) NCN2  - 
(10-15)  NCN  - 

(16-25)  JIC  - 


(26-35) SIGMA0- 


Present  crack  tip  node 
Ultimate  crack  tip  node 

Number  of  iterations  to  release  spring  (5  should  be 
adequate) 

Average  value  of  the  J-integral  paths  at  which  crack  is 
to  be  grown  (Note:  Crack  will  not  grow  unless  JIC  is 
exceeded ) 

Flow  stress  in  Tearing  Modulus  Equation  (2. 6. 3-1) 


2)  Elastic  Analysis  - 

Group  VIII  -  Nonl inear  Analysis  Data 
A.  Nonl inear  Analysis  Control 
Card  1 


If  this  card  is  left  blank  or  NINC  =  1,  the  program  will  default  to  an  APES 
run  (no  further  data  required) .  The  user  is  suggested  to  try  his  model  in 
this  mode  before  attempting  a  PAPST  run.  Cost  for  an  elastic  run  may  be  one 
to  two  orders  of  magnitude  less  expensive. 


3)  Wavefront  Reordering  -  This  is  so  the  user  can  change  element  ordering 
after  the  fact,  (i.e.,  "grid  too  big  for  progran")  without  changing 
input  model  data. 


Automatic  resequencing  of  the  elements  by  the  program  is  still  in 
developnent.  However,  the  user  may  manually  define  a  resequencing.  This  is 
particularly  useful  if  there  are  a  large  number  of  spring  elements  since 
they  are  generated  at  the  end  of  the  area  element  list.  To  use  this 
facility: 

Group  I  -  Preliminary  Data 
Card  2: 


(66-70)  IWAVE  -  option  for  wavefront  reordering 
*  0  no  change  (default) 

®  1  user  specified  sequence 


£&N3dXa  INIMNN^aoi  iv 
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Group  VII  -  Additional  Data 
D.  Wavefront  Reordering 


CARDS:  IN  (I) 

FORMAT  (1615) 

TERMINATE  WITH  A  FULLY  BLANK  CARD 

IN (I)  are  lists  of  element  numbers  that  refer  to  the  position  of  the 
element  in  the  input  stream.  Springs  are  added  on  at  the  end  of  the  area 
element  input  and  are  numbered  in  order  of  appearance  in  input.  The 
initial  value  (and  default)  of  this  list  is  reflected  in  the  output  in  the 
nodal  indices  section. 

* 

The  user  may  specify  as  many  (up  to  16)  elements  on  a  card  as  he 
wishes.  The  solution  procedure  will  be  modified  internally  to  use  the 
specified  list.  There  is  no  effect  on  the  output  sequencing. 

Unspecified  elements  in  the  list  will  be  added,  in  order,  to  the  end  of 
the  user  specified  list. 

This  solution  may  be  preferred  over  modification  of  the  element  order 
if  the  user  gets,  or  suspects  to  get,  the  error  message:  "GRID  TOO  LARGE 
FOR  PROGRAM" 

4)  Convergence  Criteria 

Group  VIII  -  Nonlinear  Analysis  Data 

Card  1  C1,C2 

The  equations  for  convergence  criteria  have  been  slightly  modified,  see 
Section  2.0.3.  As  aresult,  the  default  values  were  modified.  In  the  new 
code.  Cl  defaults  10-^  and  C2  defaults  to  10  . 


fS 
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APPENDIX  B  -  NOTES  ON  APES  AND  PAPST  CODES 


Hi  is  manual  docunents  the  current  approaches  and  algorithms  contained 
in  PAPST.  During  the  course  of  this  work  several  areas  of  concern  with 
respect  to  the  previous  versions  of  APES  and  PAPST  appeared.  A  few  coding 
errors  in  both  APES  and  PAPST  were  found.  Some  portions  of  the  codes  were 
found  to  have  questions  on  interpretation.  Many  improvements  to  the  code 
were  made  for  consistency  and  clarity.  11113  Appendix  addresses  these  areas. 

Topics 

1)  Errors  -  APES 

2)  Errors  -  PAPST 

3)  Warnings  on  Interpretation  -  Previous  Codes 

4)  Improvements  -  PAPST 

5)  Limitations  of  the  new  PAPST 


1)  Errors  -  APES 


a)  Variable  NOPT  or  ISTRN 


This  variable  changes  value  several  times  in  the  program.  The 
subroutine  BB  was  generated  from  one  section  of  the  progran  and  is  called 
from  another  where  the  designation  has  changed.  It  is  used  four  times  in 
the  subroutine.  In  order,  the  changes  are: 


N0PT.GE.2 

N0PT.EQ.2 

(x2)NOPT.NE.3 


becomes  N0PT.LE.1 
becomes  NOPT. EQ. 1 
becomes  NOPT.NE. 0 


There  may  be  other  subroutines  where  a  similar  problem  occurs,  so  the 
user  may  wish  to  review  all  references  to  this  variable. 

b)  Subroutine  FR0N1 

In  the  coding  for  inclined  or  rotated  nodes,  there  is  an  error  in  the 
sin/cos  components: 

SINA*COSA  becomes  SINA*COSB 


These  changes  have  been  implemented  in  the  versioon  of  APES  at  DTNSRDC 
2)  Errors  -  PAPST 

a)  Subroutine  B1  -  variable  interchange 

In  the  determination  of  the  enriching  functions,  the  definitions  of  two 
of  the  variables  were  interchanged.  The  code  should  read: 

TEMP8  =  TEMP8  +  DNF(2,I) *UX1 (I) 

TEMP9  =  TEMP9  +  DNF(1 , I) *UY1 (I) 

For  enriched,  large  displacement  problems,  this  error  has  an  effect  on 
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the  results 


b)  Thermal  Strains  -  DSTR 

In  plane  strain,  the  thermal  strains  depend  on  the  current  material 
property  matrix.  This  was  not  properly  handled.  See  Section  2.5.4  for  the 
correct  equations.  Since  the  material  matrix  is  formulated  differently  here 
than  in  the  latest  program,  no  simple  fix  is  described.  These  errors  have 
been  corrected  in  the  latest  version  of  PAPST  available  at  DTNSRDC. 

3)  Warnings  on  Interpretation  -  Previous  Codes 

This  section  covers  areas  where  the  user  may  have  a  misunderstanding 
about  what  the  progran  is  doing.  This  misunderstanding  may  have  an  impact 
on  how  the  user  interprets  his  results. 

a)  Subroutine  ENSTR-APES 

This  subroutine  calculates  the  crack  tip  loading  for  crack  face 
pressure.  We  have  been  unable  to  reconstruct  the  coding  in  this  section. 
However,  we  have  no  evidence  that  it  is  wrong.  The  user  may  want  to  compare 
the  crack  tip  contribution  in  the  present  subroutine  PRESS  to  that  in  ENSTR. 

b)  Jaumann  Stress  Rate  and  Enriched  Elements  -  PAPST 

The  rotation  components  used  in  the  computation  of  the  stress  rate 
correction  term  were  not  updated  for  the  enriched  contributions.  This  will 
effect  the  answers  slightly  and  are  now  included  in  the  new  program. 

c)  Hydrostatic  Expansion  -  PAPST 

We  have  found  in  .^viewing  output  for  the  previous  version  of  PAPST  that 
in  plastic  zones,  the  hydrostatic  stress  and  strain  do  not  satisfy  the 
equation: 


yy 


azz} 


(rajy 


+  e  +  e  ) 

yy  zz‘ 


(B-l) 


This  is  believed  to  be  a  numerical  inaccuracy  inherent  in  the  former  code 
and  is  satisfied  in  the  current  code. 

d)  Subroutine  OUTPUT  -  Effective  Stress  -  PAPST 

The  average  effective  stress  output  is  the  average  of  the  effective 
stresses.  It  is  not  the  effective  stress  of  the  averaged  stresses  printed 
on  the  same  line.  This  may  be  confusing  to  the  user  who  is  trying  to  use 
the  averaged  stresses  in  conjunction  with  the  effective  stress  (see  Section 
4  d) . 


>’fiNa<4Y3  iNlMNMiAon  iw  nj'innnwj^M 


80 


/k  Arthur  D.  Little,  Inc 


e)  J-Integral  Path  Values  -  PAPST 


Hie  path  taken  for  the  J-integral  is  taken  within  an  element  at  the 
edge.  The  element  is  selected  by  sorting  through  the  nodal  indices. 
Rearrangement  of  the  indices  can  result  in  a  change  in  the  elements  chosen 
for  a  given  path,  and  as  a  result,  in  a  numerically  different  J-integral 
value  (see  Section  4  f) . 

4)  Improvements  -  PAPST 

a)  Subroutine  PRESS 

This  subroutine  was  changed  from  an  edge  integral  to  a  line  integral, 
simplifying  the  code.  The  crack  tip  loading  term  has  been  generalized  to 
handle  plastic  as  well  as  elastic  singularities. 

b)  Deformation  Gradient  Contributions 

The  deformation  derivatives  (U,  y  and  V,  x)  are  now  anoothed  to  the 
nodes  in  the  same  manner  as  the  strain  components.  The  enriched  terms  have 
been  included  in  the  quadrature  point  calculations,  so  that  the  Jaumann 
Stress  Rate  is  calculated  consistently. 

c)  Yield  Surface  Modification 

The  yield  surface  is  updated  only  for  converged  increments.  The  scaling 
algorithm  for  increments  crossing  the  yield  surface  has  been  substanti tally 
revised.  This  seems  to  reduce  the  problem  of  yield  surface  "wandering" 
during  iterations.  Equation  B-l  also  seems  to  be  enforced  more  accurately. 

d)  Subroutine  CUTPT2  -  Effective  Stress 

The  effective  stress  listed  with  the  averaged  and  principle  stresses  is 
calculated  from  the  averaged  stresses.  This  keeps  the  output 
self-consistent. 

e)  Stress  Intensities 

The  stress  intensity  coding  has  been  completely  redone.  The  elastic  and 
plastic  singularities  now  use  forms  of  the  same  equations. 

f)  J-Integral,  Path  Values 

The  J-Integral  has  been  rewritten  so  that  it  uses  the  averaged  nodal 
stress  and  strain  values.  This  removes  any  sequencing  dependence  and 
appears  to  give  more  consistent  results. 

g)  Thermo-Plasticity 

The  thenmo-plastic  analysis  has  been  re-written  to  fully  incorporate  all 
of  the  plasticity  effects. 
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5)  Limitations  of  the  New  PAPST 

a)  Staoothing  Functions 

Comparison  of  nearby  quadrature  point  stress  levels  with  the  snoothed 
nodal  values  has  brought  up  a  question  of  adequacy  of  the  smoothing 
functions  in  regions  of  high  plasticity.  It  is  felt  that  some  of  the 
divergence  of  the  J-integral  path  values  may  be  a  result  of  this  limitation. 

b)  Orthotropic  Plasticity 

The  program  is  designed  to  do  orthotropic  elasticity  or  isotropic 
plasticity.  While  the  program  will  accept  orthotropic  values  for  a 
plasticity  problem'-,  the  user  is  warned  against  doing  this.  The  program  may 
or  may  not  rim  if  the  user  does  specify  these  types  of  parameters.  The 
validity  of  the  results,  if  any,  could  only  be  determined  by  the  user.  We 
have  not  tested  this  case  as  it  constitutes  an  incorrect  application  of  the 
program. 
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Appendix  C  -  JCL  and  Usaqe 


n^n  • 2 *  fl.of  PAPST  has  been  broken  into  two  separate  programs. 
SPAPST  is  the  source  code  for  the  analysis  program.  SPPOST  is  th® 
source  code  for  the  plotting  package.  The  following  files  are  of” 
interest  to  the  user: 


TAPE  5  - 

TAPES  - 
TAPE 4 7  - 

TAPE 4 8  - 


Local  name  for  the  input  file 
Local  name  for  the  output  file 

Binary  file  passed  between  PAPST  and  PPOST  with 
information  for  plotting 

Plotting  file  from  PPOST  to  be  oassed  to  plottip.q 
hardware 


5*1  gives  an  example  of  how  a  typical  job  control  file  (JCL) 
might  be  set  up  for  a  CDC  machine. 
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TABLE  C.l 


Example  JCL  Run  File 


JCL 


COMMENTS 


/JOB 

/USER 

/CHARGE 

GET  i  TAPE3=IFN . 

REWI HD t TAPES. 

GET  •>  f ' A P S T B • 

FAPST3. 

GOTO » 1  DONE . 

E/IT- 
!.  DOME  »  !  . 

• EUiNB t TAPE 6  . 
REPLACE >  TAPE6  =  OFN . 
::-.UIND>TAPE47. 

GET » PF'QSTB  » 

REWIND >TAPE48. 
REPLACE  ? TAPE48=F'FN ♦ 

GO" - 2 DONE  * 

2 I T  , 

2  DO  ME*. It* 

''!  AY F ILE.  ?  DFN  ♦ 

REPLACE  r DFN ♦ 

•  UR 


IF N “INPUT  FILE  NAME 

PAPSTB  =  BINARY  CODE  FOR  PAPST 
EXECUTE  PAPST. 


OFN=OUTPUT  FILE  NAME. 

(IF  TAPE 47  DOESN'T  EXIST  *  DROP  TO  EXIT.) 
P P OST 3= B I AN R Y  FOR  PPOST. 

EXECUTE  P POSTS . 

P F N  =  P L 0 T  FILE  NAME. 


DFN -DAYFILE  NAME. 
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